library(lavaan)
library(readxl)
library(psych)
library(openxlsx)

# Read data
data <- read_excel("/Users/qianxiaoye/Desktop/Study2 CSR/见数study2 问卷336.xlsx")

# Define all models
# Model 1: Six-factor model (baseline)
model_six <- '
  # Measurement model
  f_ECSIR =~ ECSIR1 + ECSIR2 + ECSIR3 + ECSIR4
  f_ICSR =~ ICSR1 + ICSR2 + ICSR3 + ICSR4 + ICSR5
  f_shame =~ shame1 + shame2 + shame3 + shame4 + shame5
  f_anger =~ anger1 + anger2 + anger3
  f_evoice =~ evoice1 + evoice2 + evoice3 + evoice4
  f_prosocial =~ prosocial1 + prosocial2 + prosocial3 + prosocial4
'

# Model 2: Five-factor model A (combining shame and anger)
model_five_A <- '
  f_ECSIR =~ ECSIR1 + ECSIR2 + ECSIR3 + ECSIR4
  f_ICSR =~ ICSR1 + ICSR2 + ICSR3 + ICSR4 + ICSR5
  f_emotion =~ shame1 + shame2 + shame3 + shame4 + shame5 + anger1 + anger2 + anger3
  f_evoice =~ evoice1 + evoice2 + evoice3 + evoice4
  f_prosocial =~ prosocial1 + prosocial2 + prosocial3 + prosocial4
'

# Model 3: Five-factor model B (combining evoice and prosocial)
model_five_B <- '
  f_ECSIR =~ ECSIR1 + ECSIR2 + ECSIR3 + ECSIR4
  f_ICSR =~ ICSR1 + ICSR2 + ICSR3 + ICSR4 + ICSR5
  f_shame =~ shame1 + shame2 + shame3 + shame4 + shame5
  f_anger =~ anger1 + anger2 + anger3
  f_response =~ evoice1 + evoice2 + evoice3 + evoice4 + prosocial1 + prosocial2 + prosocial3 + prosocial4
'

# Model 4: Four-factor model (combining both pairs)
model_four <- '
  f_ECSIR =~ ECSIR1 + ECSIR2 + ECSIR3 + ECSIR4
  f_ICSR =~ ICSR1 + ICSR2 + ICSR3 + ICSR4 + ICSR5
  f_emotion =~ shame1 + shame2 + shame3 + shame4 + shame5 + anger1 + anger2 + anger3
  f_response =~ evoice1 + evoice2 + evoice3 + evoice4 + prosocial1 + prosocial2 + prosocial3 + prosocial4
'

# Model 5: Three-factor model
model_three <- '
  f_CSR =~ ECSIR1 + ECSIR2 + ECSIR3 + ECSIR4 + ICSR1 + ICSR2 + ICSR3 + ICSR4 + ICSR5
  f_emotion =~ shame1 + shame2 + shame3 + shame4 + shame5 + anger1 + anger2 + anger3
  f_response =~ evoice1 + evoice2 + evoice3 + evoice4 + prosocial1 + prosocial2 + prosocial3 + prosocial4
'

# Model 6: Two-factor model
model_two <- '
  f_antecedents =~ ECSIR1 + ECSIR2 + ECSIR3 + ECSIR4 + ICSR1 + ICSR2 + ICSR3 + ICSR4 + ICSR5
  f_consequences =~ shame1 + shame2 + shame3 + shame4 + shame5 + anger1 + anger2 + anger3 + 
                    evoice1 + evoice2 + evoice3 + evoice4 + prosocial1 + prosocial2 + prosocial3 + prosocial4
'

# Model 7: One-factor model
model_one <- '
  f_general =~ ECSIR1 + ECSIR2 + ECSIR3 + ECSIR4 + 
               ICSR1 + ICSR2 + ICSR3 + ICSR4 + ICSR5 + 
               shame1 + shame2 + shame3 + shame4 + shame5 + 
               anger1 + anger2 + anger3 + 
               evoice1 + evoice2 + evoice3 + evoice4 + 
               prosocial1 + prosocial2 + prosocial3 + prosocial4
'

# Fit all models
fit_six <- cfa(model_six, data = data)
fit_five_A <- cfa(model_five_A, data = data)
fit_five_B <- cfa(model_five_B, data = data)
fit_four <- cfa(model_four, data = data)
fit_three <- cfa(model_three, data = data)
fit_two <- cfa(model_two, data = data)
fit_one <- cfa(model_one, data = data)

# Function to get comprehensive fit indices
get_indices <- function(fit) {
  chi2 <- fitMeasures(fit, "chisq")
  df <- fitMeasures(fit, "df")
  chi2_df <- chi2/df
  cfi <- fitMeasures(fit, "cfi")
  tli <- fitMeasures(fit, "tli")
  rmsea <- fitMeasures(fit, "rmsea")
  rmsea_ci_lower <- fitMeasures(fit, "rmsea.ci.lower")
  rmsea_ci_upper <- fitMeasures(fit, "rmsea.ci.upper")
  srmr <- fitMeasures(fit, "srmr")
  
  return(c(chi2, df, chi2_df, cfi, tli, rmsea, rmsea_ci_lower, rmsea_ci_upper, srmr))
}

# Create comparison table
models <- list(
  "Six-factor model (baseline)" = fit_six,
  "Five-factor model A (combining shame and anger)" = fit_five_A,
  "Five-factor model B (combining evoice and prosocial)" = fit_five_B,
  "Four-factor model" = fit_four,
  "Three-factor model (CSR, emotions, responses)" = fit_three,
  "Two-factor model (antecedents vs. consequences)" = fit_two,
  "One-factor model" = fit_one
)

# Create results dataframe
results <- data.frame(
  Model = names(models),
  do.call(rbind, lapply(models, get_indices))
)

names(results) <- c("Model", "χ2", "df", "χ2/df", "CFI", "TLI", "RMSEA", "RMSEA.LB", "RMSEA.UB", "SRMR")

# Round numeric columns to 3 decimal places
results[, 2:10] <- round(results[, 2:10], 3)

# Save detailed results
write.xlsx(results, "/Users/qianxiaoye/Desktop/Study2 CSR/cfa_detailed_results.xlsx", 
           rowNames = FALSE, colNames = TRUE)

# Perform chi-square difference tests
cat("\nChi-square Difference Tests:\n")
anova_results <- anova(fit_six, fit_five_A, fit_five_B, fit_four, fit_three, fit_two, fit_one)
print(anova_results)

# Get standardized factor loadings for baseline model
cat("\nStandardized Factor Loadings for Six-factor Model:\n")
std_loadings <- standardizedSolution(fit_six)
print(std_loadings[std_loadings$op == "=~", c("lhs", "rhs", "est.std", "pvalue")]) 